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HIGHLIGHTS 


►  A  3D  multiphase  model  is  developed  to  study  the  cold  start  process  of  PEMFC  stacks. 

►  Stacks  with  more  cells  can  reach  higher  temperature  with  better  performance. 

►  Ice  formation  in  middle  cells  is  slower. 
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To  comprehensively  understand  the  cold  start  processes  of  proton  exchange  membrane  fuel  cell  (PEMFC) 
stack  which  is  important  for  the  automotive  applications,  a  three-dimensional  multiphase  PEMFC  stack 
model  is  developed  in  this  study.  The  detailed  analysis  of  the  cold  start  processes  shows  that  for  the 
stacks  with  more  cells,  the  voltage  decreases  more  slowly  due  to  the  lower  ice  formation  rates.  The 
temperature  increases  faster  for  a  stack  with  more  cells,  and  a  higher  temperature  can  be  reached  at  the 
end  of  the  cold  start  process.  No  apparent  difference  in  voltage  exists  among  the  different  individual  cells 
in  a  stack  when  the  reactant  gases  are  evenly  supplied  to  each  cell.  The  temperature  in  the  individual  cell 
in  the  middle  of  a  stack  is  higher  and  more  evenly  distributed  than  those  on  the  sides  and  single  cells, 
due  to  weakened  cooling  effect  of  the  bi-polar  plate  (BP)  on  the  membrane  electrode  assembly  (MEA), 
and  the  ice  formation  rate  is  also  lower  in  the  middle  cell.  At  a  lower  current  density,  the  ice  in  the 
cathode  catalyst  layer  (CL)  is  formed  faster  at  the  section  close  to  the  BP,  and  it  is  close  to  the  membrane 
at  a  higher  current  density. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Proton  exchange  membrane  fuel  cell  (PEMFC)  is  one  of  the 
promising  clean  and  efficient  power  sources  for  transportation 
applications,  owing  to  its  high  power  density,  low  temperature 
operation  and  zero/low  emission  [1],  Start-up  of  PEMFC  stack  from 
sub-zero  temperatures  has  been  recognized  to  be  an  essential  issue 
before  its  successful  commercialization.  Several  countries  and 
regions  have  set  specific  targets  for  cold  start  performance.  For 
example,  in  the  United  States,  the  latest  target  was  set  by  the 
Department  of  Energy  (DOE),  which  requires  unassisted  successful 
start-up  from  -40  °C  [2],  a  much  more  difficult  requirement  than 
the  one  established  in  2005.  The  European  Union  sets  several 
general  technical  targets  for  the  years  2015-2020,  including  the 
lowest  successful  cold  start  temperature  of  -25  °C,  and  well 
maintained  proton  conductivity  at  low  temperatures  (over 
10  mS  cm-1  at  -20  °C)  [3],  The  Japanese  government  has  also 
established  a  series  of  programs,  and  in  recent  years,  the  main  plan 
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has  been  shifted  from  “strategic  development”  to  “commerciali¬ 
zation  promotion"  [4],  indicating  that  the  focus  of  research  and 
development  is  now  on  the  commercially  applied  fuel  cell  stacks. 
As  for  fuel  cell  manufacturers,  the  General  Motors  Corporation  set 
the  target  of  unassisted  start-up  from  -40  °C  as  well  [5], 

A  number  of  experimental  studies  have  been  conducted  to 
investigate  the  cold  start  performance  and  characteristics  of 
PEMFCs  [6—18],  The  focus  has  been  on  the  measurement  of  general 
performance  behaviors,  performance  degradation  and  visualization 
of  ice  formation  during  cold  start  processes.  Cold  start  models  have 
also  been  developed  for  single  cells  [19—24]  and  stacks  [25—27]  in 
an  effort  to  better  understand  the  start-up  process.  Single  cell 
analytical  models  have  been  mainly  targeted  to  discover  the  rela¬ 
tionship  between  the  design/operating  parameters  and  cold  start 
performance  [19-21],  and  most  of  these  models  have  only 
considered  either  individual  components  (e.g.  cathode  catalyst 
layer  (CL))  [19]  or  simplified  cells  (one-dimensional  models) 
[20,21  ].  Multiphase  multi-dimensional  numerical  simulations  have 
been  carried  out  for  single  PEMFCs  to  investigate  more  detailed 
transport  processes  [22-24],  Mao  et  al.  [22]  found  that  ice  in  the 
cathode  CL  appears  first  at  the  flow  channel  inlet  region  and  extends 
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Nomenclature 

desb 

desublimation 

eff 

effective 

A 

stack  geometric  area,  m2 

ele 

electronic 

c 

molar  concentration,  mol  m  3 

equil 

equilibrium 

cp 

specific  heat,  J  kg-1  K-1 

EOD 

electro-osmotic  drag 

D 

mass  diffusivity,  m2  s  1 

evp 

evaporation 

EW 

equivalent  weight  of  membrane,  1100  kg  kmol-1 

f 

frozen 

F 

Faraday’s  constant,  96,487  C  mol-1 

fl 

fluid  phase 

h 

latent  heat,J  kg-1;  heat  transfer  coefficient,  W  m  2  I<  1 

fmw 

frozen  membrane  water 

I 

current  density,  A  cm  2 

fusn 

fusion 

j 

reaction  rate,  A  m  3 

FPD 

freezing  point  depression 

k 

thermal  conductivity,  W  m  1  K-1 

g 

gas  phase 

K 

permeability,  m2 

GDL 

gas  diffusion  layer 

m 

mass  flow/transfer  rate,  kg  s-1 

h2 

hydrogen 

M 

molecular  weight,  kg  kmor1 

h2o 

water 

P 

pressure,  Pa 

i 

the  ith  components  or  the  ith  cell  in  a  stack 

Q 

heat  flux,  W  nrr2 

ice 

ice 

Q 

heat  transfer  rate,  W 

in 

inlet 

RH 

relative  humidity 

init 

initial  condition 

s 

volume  fraction 

ion 

ionic 

s 

source  terms,  entropy,  J  kmol  1  K  1 

iq 

liquid  water 

t 

time,  s 

m 

mass  (for  source  term) 

T 

temperature,  K  or  °C 

melt 

melt 

Tt 

velocity,  m  s-1 

mem 

membrane 

x 

mole  fraction 

N 

total  number  of  cells  in  a  stack,  normal  condition 

Y 

mass  fraction 

nf 

non-frozen 

nmw 

non-frozen  membrane  water 

Creek  letters 

02 

oxygen 

a 

transfer  coefficient 

out 

outlet 

£ 

porosity 

pc 

phase  change 

V 

overpotential,  V 

ref 

reference  state 

i 

interfacial  drag  coefficient 

sat 

saturation 

K 

electrical  conductivity,  S  m  1 

si 

solid  phase 

X 

water  content  in  ionomer 

stk 

stack  characteristic 

p 

dynamic  viscosity,  kg  m  1  s-1 

surr 

surroundings 

? 

stoichiometry  ratio 

T 

energy  (for  source  term) 

p 

density,  kg  m  3 

u 

momentum  (for  source  term) 

electrical  potential,  V 

vp 

water  vapor 

(0 

volume  fraction  of  ionomer  in  catalyst  layer 

wall 

surrounding  wall  of  the  single  cell  or  stacks 

0 

intrinsic  value 

Subscripts  and  superscripts 

1-i 

liquid  water  to  ice  (vice  versa) 

a 

anode 

n-f 

non-frozen  membrane  water  to  frozen  membrane 

act 

activation 

water  (vice  versa) 

BP 

bi-polar  plate 

n-i 

non-frozen  membrane  water  to  ice 

c 

cathode 

n-v 

non-frozen  membrane  water  to  vapor  (vice  versa) 

CL 

catalyst  layer 

v-i 

vapor  to  ice 

cond 

condensation 

v-1 

vapor  to  water  liquid  (vice  versa) 

toward  the  outlet  region  gradually,  the  start-up  current  density 
influences  the  water  uptake  potential  of  the  membrane,  and  a  lower 
current  density  is  beneficial  for  water  uptake.  Meng  [23]  predicted 
that  the  cold  start  process  obtains  a  benefit  from  the  higher  gas  flow 
rates  in  the  cathode  flow  channel.  Jiao  and  Li  [24]  found  that  it  is 
favorable  to  increase  the  ionomer  fraction  in  the  cathode  CL,  and 
a  thinner  membrane  enhances  the  water  uptake  rate. 

In  order  to  supply  sufficient  power  for  vehicles,  individual 
PEMFCs  are  often  assembled  in  series  to  form  stacks.  Interactions 
between  the  individual  cells  in  a  stack  become  significant,  making 
the  cold  start  characteristics  of  PEMFC  stacks  different  from  single 
cells  [16],  Sundaresan  and  Moore  [25]  developed  a  layered  cold 
start  stack  model,  which  considered  the  water  phase  change  and 
the  thermal  effects,  and  this  model  predicted  the  temperature  of 
different  cells  in  a  stack.  A  one-dimensional  model  developed  by 


Khandelwal  et  al.  [26]  also  obtained  the  temperature  distribution  in 
a  PEMFC  stack.  Based  on  this  model,  several  assisted  start-up 
strategies  were  tested  and  compared.  Ahluwalia  and  Wang  [27] 
developed  a  two-dimensional  cold  start  stack  model,  and  the  ice 
formation  was  considered  in  this  model.  It  was  found  that  a  high 
current  density  is  favorable  for  rapid  cold  start  of  PEMFC  stacks,  and 
the  metal  bi-polar  plates  (BPs)  are  better  than  graphite  to  improve 
the  cold  start  abilities  due  to  the  lower  heat  capacity  of  metal.  The 
PEMFC  cold  start  related  studies  were  also  reviewed  recently  in 
detail  by  Meng  and  Ruan  [28]  and  Jiao  and  Li  [29], 

As  mentioned  above,  most  of  the  previous  cold  start  models 
only  considered  single  cells  and  one-dimensional/two-dimensional 
approach  to  stack  level.  To  comprehensively  understand  the  PEMFC 
stack  cold  start  processes,  a  three-dimensional  multiphase  PEMFC 
stack  model  is  needed.  In  this  study,  a  three-dimensional 
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multiphase  PEMFC  stack  model  is  developed,  and  the  cold  start 
processes  of  PEMFC  stacks  are  investigated  numerically.  This  article 
is  organized  as  follows:  the  stack  model  is  presented  first,  followed 
by  the  results  and  discussion,  and  finally  the  conclusion  is  given. 

2.  Stack  model 

2.1.  Physical  problem 

The  computational  domain  of  this  model  is  based  on  practical 
PEMFC  stacks  with  typical  components,  including  the  BPs  with  flow 
channels  machined  on  them  and  the  membrane  electrode  assem¬ 
blies  (MEAs).  To  save  the  computational  time,  a  50  mm  straight 
segment  of  each  cell  in  the  stacks  are  considered,  which  consists  of 
one  of  the  parallel  channels  as  well  as  the  certain  part  of  BPs,  gas 
diffusion  layers  (GDLs),  CLs  and  membranes.  The  stack  component 
properties  and  operating  conditions  are  shown  in  Table  1,  and  the 
computational  domain  and  mesh  for  a  three-cell  stack  considered  in 
this  study  are  shown  in  Fig.  1.  With  more  computational  power  and 
time,  stacks  with  more  individual  cells  and  larger  reaction  areas  can 
also  be  simulated  by  using  this  model  with  general  formulations. 

2.2.  Assumptions 

A  number  of  assumptions  are  made  for  this  stack  model: 

(1)  The  flow  is  laminar  due  to  the  low  Reynolds  number. 

(2)  Isotropic  material  properties  are  considered. 

(3)  All  the  gas  species  are  ideal  gas. 

(4)  Water  exists  in  the  ionomer  as  non-frozen  membrane  water 
and  frozen  membrane  water,  in  the  pores  of  CL  and  GDL  as 
vapor,  liquid  and  ice,  and  the  amount  of  liquid  water  in  flow 
channel  is  fixed  at  zero  because  it  can  be  removed  quickly  in 
the  short  and  straight  flow  channels  considered  in  this  study. 

(5)  The  product  water  from  the  electrochemical  reaction  is 
absorbed  by  the  ionomer  in  the  cathode  CL  immediately,  thus 
in  the  non-frozen  membrane  water  state. 

(6)  The  different  phases  (gas,  liquid  and  solid)  are  at  local  ther¬ 
modynamic  equilibrium  such  that  they  share  the  same 
temperature  at  the  same  grid  point,  therefore  a  single  energy 
conservation  equation  is  solved  in  this  model. 

2.3.  Conservation  equations 

Mass  conservation  of  gas  mixture: 

K<i -*-*.)*) +v-(**.)-s.  (,) 

where  Sm  (kg  m  3  s-1)  is  the  source  term: 


Table  1 

Stack  properties  and  operating  conditions. 


Parameter 

Value 

Channel  length;  width;  depth; 

50;  1.0;  1.0;  1.0  mm 

rib  width 

Thicknesses  of  membrane 

0.051;  0.01;  0.2  mm 

(Nafion  112);  CL;  GDL 

n _ _  ™  on  =  1 980:  1 000: 1 000: 

Densities  of  membrane:  CL:  ’  ’ 

GDL;  BP 

1000  kg  m-J 

Specific  heat  capacities  of 

(cp)mem,CL,GDL,BP  =  833;  3300;  568; 

membrane;  CL;  GDL;  BP 

1580J  kg-1  K"1 

Thermal  conductivities  of 

kmemCLGDLBP  -  0.95;  1 .0;  1 .0; 

membrane;  CL;  GDL;  BP 

20  Wm-1  K"1 

Electric  conductivities 

Kcl,gdl,bp  =  300;  300;  20, 000  S  m”1 

of  CL;  GDL;  BP 

Volume  fraction  of  ionomer  (w) 

0)  =  0.2;  eclgdl  =  0.3;  0.6 

in  CL,  and  the  porosities  (e) 
of  CL;  GDL 

Intrinsic  permeabilities  of  CL; 

Kq.  =  6.2  x  10"13  m2; 

GDL 

XgDL  =  6.2  x  1 Q-  '2  m2 

Current  densities 

I  =  0.05  A  cm-2  or  0.1  A  cm"2 

(galvanostatic) 

Stoichiometry  ratio 

Uc  =  2.0 

Relative  humidities  of 

rhSv  =  o 

inlet  gases 

Inlet  gas  and  surrounding 

7™.  =  rsUrr  =  -20  °C 

temperatures 

Pressure  at  outlets 

p°f  =  1  atm 

Initial  stack  temperature 

Trtkjnit  =  Tsurr  =  -20  °C 

Heat  transfer  coefficients 

h  =  50  W  m"2  K"1 

Initial  ice  volume  fraction 

0 

Initial  non-frozen  water 

5;  0 

content  in  the  membrane 
and  CLs  and  frozen  water 
content  in  the  membrane 

oxygen  consumption  in  cathode  CL;  and  the  third  term  accounts  for 
the  phase  change  related  to  water  vapor: 


Sh2 


So2 


I 

I 


-|^Mh2  (in  anode  CL) 

0  (in  other  zones) 

-Jm02  (in  cathode  CL) 
0  (in  other  zones) 


Svp 


-Sv_,  -  Sv-i  +  Sn-vMH2o  (in  CL) 
— Sv_i  -  Sy-i  (in  GDL) 


(3) 

(4) 

(5) 


Sm  =  Sh2+S02+SvP  (2) 

The  first  term  on  the  right  hand  side  of  Equation  (2)  represents  the 
hydrogen  consumption  in  anode  CL;  the  second  term  is  for  the 


The  corresponding  source  terms  for  water  condensation/evapo¬ 
ration  (Sv-i,  kg  m  3  s-1),  desublimation  (Sv-i,  kg  m  3  s-1), 
and  ionomer  hydration/dehydration  due  to  water  vapor  (Sn.v, 
kg  m-3  s-1)  are 


ic(l  -Siq-Sj 


(pgXVp  -  Psat  j  Mh2o 


(if  PgXvp  >  Psat) 


(pgXvp  -  Psat)MH2o 


(if  PgXVp  <  Psat  j 


(6) 


(if  T  <  TN  +  TfpD) 
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V 


{Tdesb  £  ( 1  ^lq  ^ice) 

0  (if  PgXvp  <  psat) 


(pgXVp  -  Psat^Mn2o 

«r 


(if  PgXVp  >  pSaL  j 


(if  t  <  rN  +  rFPD) 


(ifT>TN  +  TFPD) 


(7) 


Momentum  conservation  of  gas  mixture: 


9  f 

9t  Vt(l-slq-sice 

-  Vpg  +  PgV-( v( 


pgtig'tig 

E2(l  -Siq-Si( 


<i-slq-sice)r  v<i-slq-sice); 


Dfff  =  V5(l-slq-sice)15  (12) 

Mass  conservation  of  liquid  water: 

a(gS|qPlq)  +  v  =  v.  (PlqD|qVSlq)  +  Slq  (13) 

where  the  source  term  Siq  (kg  m  3  s-1)  includes  the  water  evapo¬ 
ration/condensation  and  freezing/melting: 

V  =  Vi  ~  Sh  (14) 

where  Sh  (kg  m-3  s-1)  represents  water  freezing/melting: 


and  the  source  term  in  this  equation  is: 

!-^t?g  (in  CL  and  GDL) 

Ks  '  (10) 

0  (in  other  zones) 

where  the  pg  (kg  m  1  s  ’)  is  the  dynamic  viscosity  of  the  gas 
mixture,  7?g  (m  s-1)  the  gas  phase  superficial  velocity  vector,  and  I(s 
(m2)  the  gas  phase  permeability. 

Mass  conservation  of  gas  species  (i:  hydrogen,  oxygen  and  water 
vapor): 


7fusn£SlqPlq  (>f  T  <  TN  +  TFPd) 
-Tmelt^icePice  (if  T  >  TN  +  TFPD) 


(15) 


The  water  phase  change  is  also  related  to  the  Gibbs-Thompson 
undercooling  effect,  and  more  detailed  discussions  on  the  water 
phase  change  considered  in  this  model  can  be  found  in  Ref.  [24], 
Mass  conservation  of  ice: 


9(£SicePice)  =  ^ 
dt  " 


(16) 


I  «1  -  Slq  -  Sice) Pg Vi)  +  V-  (pgUgn)  =  V-  (pg DfvY{)  +  S, 

(11) 

and  the  source  term  in  the  above  equation  represents  the  source 
terms  of  hydrogen,  oxygen  and  water  vapor  conservation  equa¬ 
tions,  as  given  in  Equations  (3)— (5).  The  effect  of  ice  blockage  on  the 
diffusion  of  reactant  in  electrode  is  also  considered: 


where  the  source  term  considering  the  various  water  phase 
changes  is 


Sice 


Sv_i  +  S]_j  +  Sn_jMH2o  (in  CL) 
Sv-i  +  Sw  (in  GDL) 


(17) 


where  Sn-i  (kg  m  3  s  ')  represents  the  process  that  non-frozen 
membrane  water  in  CL  freezes  to  ice: 


Fig.  1.  Computational  domain  and  mesh  (individual  cells  in  the  three-cell  stack  are  marked  as  Cell  1,  Cell  2  and  Cell  3  from  the  cathode  side  to  the  anode  side). 
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(4f  -  *sat)  (l  -  Slq  -  Sice)  (if ^  ^ 


Normal  Operating  Condition: 

In  Membrane 


Non-frozen 
I  Membrane  Water  I 


This  is  based  on  the  assumption  that  in  CL,  the  non-frozen 
membrane  water  may  freeze,  and  the  formed  ice  appears  in  the 
pores  of  CL  directly  [24], 

Mass  conservation  of  non-frozen  membrane  water: 
Pmem9(^nf)  _  faemy,  /t,5n .  \  ,  c  Mg') 

EW  m  EW  VW  nmwV/nf  j  +  Jnmw 

where  the  source  term  considering  the  water  production  from  the 
electrochemical  reactions,  electro-osmotic  drag  (EOD)  effect  and 
phase  changes  is 

,_f  (in  membrane) 

Snmw  =  l  ll  -  Sn-v  -  SnA  +  SE0D  (in  cathode  CL)  (20) 
l  -Sn-v  -  Sn_;  +  Shod  (in  anode  CL) 

Mass  conservation  of  frozen  membrane  water: 


J(wAf) 


-  =  Sfav 


Q  =*vpPg  +  2 


(=  4.837  (if  T<  223.15  K) 

i=  [- 1.304 +0.01479T- 3.594  x  10-5T2]-1  (if  223.15  K  <  T  <  Tn) 
>Anf  (if  T  >  Tn) 


Sfmw  =  Sn.f 

(22) 

(23) 

seod  - 

(24) 

In  membrane,  the  phase  change  between  frozen  and  non-frozen 
water  is  considered  (Equation  (23)).  Seod  (kmol  m-3  s_1)  represents 
the  water  transport  due  to  the  electro-osmotic  drag  effect.  The 
schematic  of  the  water  phase  change  considered  in  this  model  is 
also  shown  in  Fig.  2  [24], 

With  the  membrane  hydration/dehydration  considered  in 
Equations  (8)  and  (18),  the  non-frozen  and  frozen  water  contents 
(Anf.f),  equilibrium  water  content  (Aequil),  water  activity  (a)  and 
saturation  water  content  (ASat)  are  defined  as 

.  EW 

/lnf.f  =  ~ - cnf.f 

(25) 

,  f 0. 043  +  17. 81a-39.85a2+36.0a3  (ifO 

et>uil  \  14.0  +  1.4(0-1)  (ifl  <a<3) 

<a<!)  (26) 

The  ionic  conductivity  is  therefore  a  function  of  non-frozen 
water  content  rather  than  the  total  water  content  (the  effect  of 
water  freezing  on  membrane  conductivity  is  therefore  considered): 

,qon  =  (0.5139Anf- 0.326) exp  [1268(3^-^]  (29) 

Electron  transport: 

0  =  V-(icfeV0ele)+Sele  (30) 

Ion  transport: 


0  =  V'(KfonV0ion)  +  Sjon 

The  source  terms  in  the  above  two 
reaction  rate  in  anode  and  cathode  CLs: 


(31) 

equations  depend  on  the 


(  — ja  (in  anode  CL) 
Seie  =  ^  jc  (in  cathode  CL) 
(in  other  zones) 


■If 

■If 


(in  anode  CL) 

■Sion  =  {  -Jc  (in  cathode  CL) 
(in  other  zones) 

where  the  reaction  rates  are  defin 
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as  an  example). 


A  linear  relationship  between  the  ice  formation  and  reaction  rate  is 
considered  [24],  therefore,  if  the  CL  is  fully  blocked  by  ice,  the 
electrochemical  reaction  is  stopped  (the  cold  start  process  is  failed). 
Energy  conservation  equation: 

It  {(pCP)ZT)  +  v'  ((pCP)f  -  V-  (kfl®vr)  +  Sr  (36) 

where  the  last  term,  Sj  (W  m-3),  is  the  source  term: 


The  heat  generation  includes  the  ohmic  heat  generated  from  the 
electric  current  though  the  stacks,  activational  heat  generated  from 
the  electrochemical  reaction,  reversible  heat  due  to  the  entropy 
change,  and  the  latent  heat  due  to  the  water  phase  change.  The 
latent  heat  (Equation  (38))  includes  the  corresponding  heat 
generation  and  consumption  caused  by  the  water  phase  change 
considered  in  this  study  [30],  More  details  of  the  model  formulation 
can  also  be  found  in  Refs.  [24,31], 

2.4.  Boundary  and  initial  conditions 

For  single  cell  models,  the  voltage  and/or  current  density  are 
often  set  on  the  two  end  surfaces  of  the  end  plates  (e.g.  [31,32]). 
However,  with  more  than  one  cells  between  the  end  plates,  the 
electrical  boundary  condition  is  needed  for  each  single  cell. 

Therefore,  as  shown  in  Fig.  3,  internal  electrical  boundary 
condition  is  specified  in  the  middle  of  each  BP.  Two  virtual  surfaces 
are  assumed  existing  at  the  middle  of  each  BP,  and  the  two  surfaces 
belong  to  the  separated  cells  (one  for  the  end  of  the  anode  of  one 
cell,  and  another  for  the  end  of  the  cathode  of  the  next  cell). 


'jahactl  +  llV^elell2^  +  ll^ionl^fon  +  V  (in  anode  CL) 

+Jcbactl  +  IWelell^de  +  IW.onll^fon  +  Spc  (in  Cathode  CL) 
l|V</>eiel|2^+Spc  (in  GDL) 

IWelell2^  (inBP) 

II' V0iOnl|2<on  +  V  (in  membrane) 

,  0  (in  other  zones) 


(37) 


where  the  latent  heat,  Spc  (W  m  3),  is  defined  as 

IfrcondVi  +  (ftcond  +  Ksn)Sv-i  +  hfxlsnSM  (in  GDL) 
hcond(Vl  -  Sn-vM^o)  +  (^cond  +  hfusnlVi  +  hfusn(Sl-i  +  Sn-iMH2o)  (in  CL) 
hfusnSn-fMthO  (in  membrane) 

0  (in  other  zones) 


(38) 
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Fig.  4.  Comparison  between  the  model  predictions  and  experimental  data  (a:  from 
-10  °C  at  0.07  A  cm2;  b:  from  -5  °C  at  0.3  V)  [37],  The  stoichiometry  ratio  is  2.0,  the 
initial  water  content  is  6.2,  the  initial  ice  volume  fraction  is  0,  and  the  inlet  relative 
humidity  is  0  [37], 


Considering  the  fact  that  each  single  cell  works  at  the  same  current 
density  (connected  in  series),  at  the  ends  of  the  anodes  of  the  single 
cells  (on  the  virtual  surfaces  for  anode  and  the  anode  end  plate 
surface),  the  current  density  is  defined  to  be  the  same  for  each  cell 
to  simulate  the  galvanostatic  (constant  current  density)  start-up 
condition: 

/ f  =  j|  =  ...  =  ]j>  =  •••  -  I*  =  Constant  (39) 

where  l  (A  cm-2)  is  the  current  density;  the  subscripts  1,  2  and  i 
represent  Cell  1,  Cell  2  and  Cell  i  (the  ith  cell)  in  a  stack,  respec¬ 
tively;  N  represents  the  total  number  of  cells;  and  the  superscript 
a  represents  anode.  The  electronic  potential  is  defined  at  the  ends 
of  the  cathodes  of  the  single  cells  (on  the  virtual  surfaces  for 
cathode  and  the  cathode  end  plate  surface): 


Time,  s 


Fig.  5.  Evolution  of  average  voltages  and  volume  averaged  ice  volume  fraction  in  the 
cathode  catalyst  layers  for  single  cell,  two-cell  stack  and  three-cell  stack  at  the  current 
densities  of  (a)  0.1  A  cm~2  and  (b)  0.05  A  cirr2.  The  initial,  ambient  and  inlet  gas 
temperatures  are  all  -20  °C,  the  initial  water  content  is  5,  and  the  initial  ice  volume 
fraction  is  0. 


‘/’ele.l  =  0ele,2  =  -  =  0ele,i  =  •  =  <Pele:N  =  0  (40) 

where  </>  (V)  represents  the  electrical  potential,  and  the  subscript  ele 
represents  electrons.  This  boundary  condition  was  also  used  in 
a  previous  stack  model  [33].  The  assumption  in  this  kind  of 


Case  no.  Case  1  Case  2  Case  3  Case  4  Case  5  Case  6 

Number  of  cells  1  1  2  2  3  3 

Current  density  (A  cm-2)  0.05  0.1  0.05  0.1  0.05  0.1 


Fig.  6.  Evolution  of  volume  averaged  stack  and  single  cell  temperatures  for  different 
current  densities  (0.1  A  cirr2  and  0.05  A  cm-2)  and  numbers  of  cell  (one  cell,  two  cells 
and  three  cells).  The  initial,  ambient  and  inlet  gas  temperatures  are  all  -20  °C,  the 
initial  water  content  is  5,  and  the  initial  ice  volume  fraction  is  0. 
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boundary  condition  is  that  the  current  density  and  voltage  are 
evenly  distributed  at  the  interfaces  between  the  individual  cells.  In 
this  study,  the  active  reaction  area  is  only  1  cm2  (single  straight 
flow  channel),  much  smaller  than  the  commercial  PEMFC  stacks  on 
the  level  of  hundreds  cm2  [34,35],  and  the  variation  in  1  cm2  was 
found  to  be  insignificant  [35],  Therefore,  it  is  reasonable  to  assume 
that  the  current  density  and  voltage  are  evenly  distributed  at  the 
interfaces. 

There  should  be  no  heat  source  or  contact  heat  transfer  resis¬ 
tance  between  the  pairs  of  virtual  surfaces  defined  above,  there¬ 
fore,  for  each  pair  of  virtual  surfaces,  the  internal  thermal  boundary 
condition  is  defined  as 


Tf  =  (41) 

<?f  =  <?i,l  (42) 

where  T  (K)  stands  for  temperature,  the  equal  temperature  repre¬ 
sents  that  there  is  no  contact  heat  transfer  resistance;  and  q 
(W  m-2)  is  the  heat  flux,  representing  the  energy  conservation 
between  the  two  virtual  surfaces. 

The  inlet  boundary  conditions  of  the  flow  channels  are  specified 
the  same  for  each  inlet  port  of  the  single  cell  in  a  stack,  based  on  the 
common  requirement  of  uniform  reactant  supply  [36].  The  mass 


»  +  * 


°c 
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)n  in  the  middle  cross  section  in  x-y  planes  (25  mm  iron 
water  content  is  5,  and  the  initial  ice  volume  fraction  is 


:  gas 


Fig.  7.  Transient  temperature  distributic 
temperatures  are  all  -20  °C,  the  initial 


■  inlets)  of  the  three-cell  stack  at  0.1 
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flow  rates  at  the  inlets  are  defined 
parameters: 


_  Pg? aW 

2Fch 


Pg^chefA 

4Fc02 


function  of  several 


(43) 


where  /ref  (A  cm-2)  is  the  reference  current  density;  (kg  m  3) 
and  Pg  (kg  m-3)  are  the  densities  of  the  gas  mixtures  supplied  to 
flow  channels;  A  (m2)  is  the  flat  area  of  reaction  surface;  and  c 
(kmol  m-3)  represents  the  reactant  concentration.  The  inlet  gas 
temperature  is  set  to  be  equal  to  the  ambient  temperature  of 
-20  °C,  the  stoichiometric  ratio  is  2.0,  and  the  inlet  relative 
humidity  is  0. 

The  sidewalls  in  this  straight  channel  stack  model  are  the 
middle  surfaces  between  the  parallel  flow  channels,  on  these  walls 
the  surface  condition  is  set  to  be  adiabatic  due  to  the  insignificant 
temperature  difference  between  the  neighboring  channels  as  pre¬ 
sented  in  the  previous  experimental  work  [14],  On  the  end  walls  of 
the  end  plates  (as  shown  in  Fig.  1 ),  a  heat  transfer  coefficient  and 
the  surrounding  temperature  are  used  to  define  thermal  boundary 
condition: 


Q  =  Mwall(Tsurr  -  Twall)  (44) 

where  Q.  (W)  represents  the  heat  transfer  rate  between  the  end 
walls  of  the  stack  and  the  ambient;  h  (W  m  2  K-1)  is  the  heat 
transfer  coefficient;  Awaii  (m2)  represents  the  surface  area  of  the  end 
wall;  Tsurr  (K)  is  the  surrounding/ambient  temperature;  and  Twan 
(K)  is  the  surface  temperature  of  the  end  wall  of  the  stack. 

For  the  initial  conditions,  the  initial  ice  volume  fraction  is  set  to 
be  0,  representing  an  effective  purging  process  before  the  next  cold 
start  cycle;  an  initial  water  content  of  5.0  is  specified  in  the 
membranes  and  CLs;  and  the  initial  stack  temperature  is  set  to  be 
-20  °C,  which  equals  to  the  ambient  temperature. 


2.5.  Numerical  procedures 

This  PEMFC  stack  cold  start  model  is  implemented  with  the 
user-defined  function  (UDF)  of  the  computational  fluid  dynamics 
(CFD)  program  FLUENT.  Finite  volume  method  is  adopted  to  dis¬ 
cretize  and  solve  the  conservation  equations.  Grid  independent 
study  has  been  carried  out  to  ensure  the  rationality  of  the  simu¬ 
lations,  and  the  whole  computational  domain  has  228,000  cells  for 
the  three-cell  stack  shown  in  Fig.  1.  The  second-order  upwind 
discretization  method  is  used.  The  pressure-implicit  with  splitting 
of  operators  (P1SO)  algorithm  is  adopted  which  is  based  on  the 
higher  degree  of  the  approximate  relation  between  the  corrections 
for  pressure  and  velocity,  and  it  is  effective  in  decreasing  the  iter¬ 
ation  times  required  for  convergence  when  solving  a  transient 
problem.  The  convergence  is  accelerated  by  an  algebraic  multigrid 
(AMG)  method  with  the  Gauss— Seidel  smoother.  The  minimum 
and  maximum  time  step  sizes  are  0.005  s  and  0.1  s  for  the  simu¬ 
lations.  Convergence  criteria  with  a  residual  of  10~8  are  used  for  all 
the  variables. 

3.  Results  and  discussion 

This  model  has  been  comprehensively  compared  with  experi¬ 
mental  results  with  reasonable  agreements  [37],  As  shown  in  Fig.  4, 
both  the  failed  (from  -10  °C  at  0.07  A  cm-2)  and  successful  (from 
-5  °C  at  0.3  V)  cold  start  processes  are  compared,  and  the  operating 
conditions  in  the  numerical  simulations  are  defined  corresponding 
to  the  experiments  [37],  The  stoichiometry  ratio  is  2.0,  the  initial 
water  content  is  6.2,  the  initial  ice  volume  fraction  is  0  in  all 
components,  and  the  inlet  relative  humidity  is  0. 


Fig.  8.  Various  heat  generation  and  loss  rates  during  the  cold  start  process  of  the  single 
cell  at  0.1  A  cm-2;  (a)  transient  variation;  (b)  percentage  distribution  for  the  complete 
cold  start  process.  The  initial,  ambient  and  inlet  gas  temperatures  are  all  -20  °C,  the 


Fig.  9.  Various  heat  generation  and  loss  rates  during  the  cold  start  process  of  the  two¬ 
cell  stack  at  0.1  A  cm-2;  (a)  transient  variation;  (b)  percentage  distribution  for  the 
complete  cold  start  process.  The  initial,  ambient  and  inlet  gas  temperatures  are 
all  -20  °C,  the  initial  water  content  is  5,  and  the  initial  ice  volume  fraction  is  0. 
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Fig.  10.  Various  heat  generation  and  loss  rates  during  the  cold  start  process  of  the 
three-cell  stack  at  0.1  A  cm-2;  (a)  transient  variation;  (b)  percentage  distribution  for 
the  complete  cold  start  process.  The  initial,  ambient  and  inlet  gas  temperatures  are 


The  results  of  this  study  are  discussed  in  four  parts:  the  first  part 
compares  the  single  cell  and  the  stacks  with  different  numbers  of 
cells;  the  second  part  compares  the  different  cells  in  the  three-cell 
stack;  the  third  part  compares  the  single  cell  to  the  middle  cell 
inside  the  three-cell  stack;  and  the  last  part  presents  the  effect  of 
current  density.  The  results  for  the  effect  of  cell  numbers  on  start¬ 
up  time  in  Ref.  [26]  showed  that  the  start-up  time  changes  from 
single  cell  to  stack,  and  when  the  cell  number  reaches  three,  the 
performance  becomes  similar  with  more  cell  numbers.  Therefore, 
only  single  cell,  two-cell  stack  and  three-cell  stack  are  considered  in 
this  study  to  save  the  computational  time,  and  the  computational 
domain  for  the  three-cell  stack  with  straight  flow  channels  is 


Fig.  12.  Evolution  of  volume  averaged  cell  temperatures  and  voltages  of  the  individual 
cells  in  the  three-cell  stack  at  the  current  densities  of  (a)  0.1  A  cm~2  and  (b) 
0.05  A  cm-2.  The  initial,  ambient  and  inlet  gas  temperatures  are  all  -20  “C,  the  initial 


shown  in  Fig.  1.  Six  cases  are  considered  in  this  study  as  listed  in 
Table  2.  All  the  simulations  are  stopped  when  the  volume  averaged 
ice  volume  fraction  in  cathode  CL  approaches  1,  because  this  is  the 
indication  of  a  completely  failed  cold  start  process.  Galvanostatic 
(constant  current  density)  start-up  mode  is  adopted  for  all  the 
cases;  the  number  of  cells  varies  from  one  to  three;  and  two 
different  current  densities  of  0.05  A  cm-2  and  0.1  A  cm  2  are 
considered.  Detailed  design  and  operating  parameters  are  given  in 
Table  1. 
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Fig.  11.  Water  content  of  diverse  states  in  different  stacks  and  single  cell  at  0.1  A  cnr2. 
The  initial,  ambient  and  inlet  gas  temperatures  are  all  -20  °C,  the  initial  water  content 
is  5,  and  the  initial  ice  volume  fraction  is  0. 
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Fig.  13.  Evolution  of  cell  voltages  and  volume  averaged  ice  volume  fractions  in  cathode 
catalyst  layers  for  the  single  cell  and  Cell  2  of  the  three-cell  stack  at  different  current 
densities.  The  initial,  ambient  and  inlet  gas  temperatures  are  all  -20  °C,  the  initial 
water  content  is  5,  and  the  initial  ice  volume  fraction  is  0. 
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3.1.  Comparison  of  single  cell  and  stacks 

Fig.  5  shows  the  comparison  of  the  average  voltage  of  each  cell 
and  volume  averaged  ice  volume  fraction  in  the  cathode  CLs  in  the 
stacks  and  the  single  cell.  It  can  be  noticed  that  the  cell  voltages  are 
almost  the  same  for  the  stacks  and  single  cell  at  the  beginning,  as 
the  ice  in  the  CLs  accumulates,  the  voltage  decreases  gradually,  and 
the  difference  among  the  stacks  and  single  cell  becomes  more  and 
more  notable.  Furthermore,  the  voltage  decreases  more  slowly  if 
the  stack  has  more  cells.  The  evolution  of  the  ice  volume  fraction  is 


also  similar  among  the  different  stacks  and  single  cell,  and  as  the 
number  of  cells  increases,  the  ice  volume  fraction  increases  more 
slowly.  At  the  lower  current  density  of  0.05  A  cm~2,  since  the  ice 
formation  rate  is  lower  than  that  at  the  higher  current  density  of 
0.1  A  cm-2,  the  difference  in  start-up  time  is  more  significant. 

The  difference  in  the  voltages  and  ice  formation  rates  are  related 
to  the  different  temperature  evolution,  as  shown  in  Fig.  6.  It  shows 
that  the  volume  averaged  temperature  increases  faster  for  the 
stacks  with  more  cells,  which  can  reach  higher  temperature  at  the 
end  of  the  cold  start  processes.  It  can  also  be  noticed  that  a  higher 
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Fig.  14.  Distribution  of  ice  volume  fraction  during  the  start-up  processes  in  the  cathode  catalyst  layers  in  x-y  planes  in  the  center  of  the  stack  (25  mm  from  the  inlets)  at  a  current 
density  of  0.1  A  cm2;  (a)  the  cathode  catalyst  layer  of  the  single  cell;  (b)  the  cathode  catalyst  layer  of  Cell  2  in  the  three-cell  stack.  The  initial,  ambient  and  inlet  gas  temperatures  are 
all  -20  °C,  the  initial  water  content  is  5,  and  the  initial  ice  volume  fraction  is  0. 


110 


Y.  Luo  et  al.  /  Journal  of  Power  Sources  224  (2013)  99-1 


current  density  is  favorable  for  increasing  the  temperature.  At  the 
current  density  of  0.1  A  cm-2,  the  single  cell  can  reach  the 
temperature  of  -16.8  °C,  while  the  two-cell  stack  reaches  -16.5  °C, 
and  the  three-cell  stack  sustained  a  bit  longer  which  reaches 
-16.2  °C.  With  a  lower  current  density,  the  temperature  difference 
between  different  stacks  and  single  cell  is  larger  due  to  the  longer 
start-up  processes.  The  cell  voltage  depends  on  the  reversible 
voltage  and  overpotential,  and  the  reversible  voltage  decreases  as 
the  temperature  increases,  while  the  overpotential  is  lower  at 
higher  temperature  and  higher  with  more  ice  (less  reaction  area). 
Obviously,  the  overpotential  plays  the  major  role  in  this  case.  In 
addition,  the  ice  formation  rate  is  lower  at  higher  temperature 
because  the  ionomer  can  absorb  more  product  water. 

The  transient  temperature  distribution  of  the  three-cell  stack  is 
shown  in  Fig.  7.  The  temperature  distribution  of  the  single  cell  and 
two-cell  stack  in  the  same  cross  section  is  not  shown  here,  the 
single  cell  condition  can  be  found  in  the  previous  work  [24],  and  the 
two-cell  stack  condition  is  similar  to  the  three-cell  stack.  During 
a  cold  start  process,  the  temperature  is  mainly  influenced  by  four 
different  kinds  of  heating  sources  as  previously  mentioned.  As 
shown  in  Fig.  7,  for  each  individual  cell  in  the  stack,  the  temperature 
increases  through  the  cell  from  the  MEA  to  the  BPs,  because  the 
main  heat  sources  are  the  ohmic  heat  from  the  membranes  and  CLs 
as  well  as  the  activational  heat  from  the  cathode  CLs.  The  temper¬ 
ature  of  the  cell  in  the  middle  is  higher  than  those  outside  because 
the  middle  cell  is  further  away  from  the  end  plate  surfaces  which 
dissipate  heat  to  the  surroundings.  The  temperature  of  the  cathode 
end  plate  is  lower  than  the  anode,  because  more  ohmic  heat  is 
generated  on  the  anode  side.  The  ohmic  heat  is  higher  on  the  anode 
side  because  the  ionomer  close  to  the  anode  is  dryer  than  the 
cathode  due  to  the  electro-osmotic  drag  effect,  resulting  in  higher 
ohmic  resistance  of  the  ionomer  close  to  the  anode. 

To  further  explain  the  temperature  distribution  in  Fig.  7,  the 
evolution  of  various  heat  generation  and  loss  rates  for  the  single 
cell  and  stacks  is  shown  in  Figs.  8—10.  Fig.  8(a)  shows  the  transient 
variation  of  the  total  heat  generation  rate  and  heat  loss  rates 
through  different  ways  for  the  single  cell,  and  Fig.  8(b)  shows  the 
corresponding  percentages  of  total  heat  generation  and  losses  for 
the  complete  cold  start  process  (0—18  s);  and  the  corresponding 
results  for  the  two-cell  stack  and  three-cell  stack  are  shown  in 
Figs.  9  and  10.  In  Figs.  8(a),  9(a)  and  10(a),  the  total  heat  generation 
rate  is  the  sum  of  the  four  heat  sources  including  the  ohmic  heat, 
activational  heat,  reversible  heat  and  latent  heat.  For  the  single  cell 
and  stacks,  the  total  heat  generation  rate  increases  slightly  during 
the  cold  start  process  as  a  combined  effect  of  the  dryness  of  the 
anode  CL  and  an  increase  of  overpotential.  As  three  ways  of  heat 
transfer  are  considered  which  include  the  heat  absorbed  by  PEMFC, 
loss  from  walls  (end  plate  surfaces)  and  loss  from  reactant  gases. 
The  tendencies  of  transient  variation  for  the  heat  transfer  are  also 
similar  among  the  single  cell  and  stacks.  As  shown  in  Fig.  8(a),  the 
proportion  of  heat  which  dissipates  to  environment  increases 
during  the  start-up  process,  correspondingly,  the  proportion  of 
heat  to  raise  the  cell  temperature  decreases.  The  increment  of  heat 
loss  rate  is  a  result  of  the  increased  cell  temperature.  Moreover,  due 
to  the  low  heat  capacity  of  hydrogen  and  air  comparing  with  the 
solid  materials,  the  heat  loss  from  the  gases  is  insignificant.  By 
comparing  the  different  stacks  and  single  cell,  it  can  be  found  that 
the  total  heat  generation  rate  increases  almost  proportionally  with 
the  cell  numbers,  and  together  with  the  fact  that  the  increased  rates 
of  heat  loss  from  the  walls  are  almost  the  same,  the  proportion  of 
heat  loss  from  the  walls  decreases  as  the  cell  number  increases 
while  the  proportion  of  heat  absorbed  by  PEMFC  increases,  as 
shown  in  Figs.  8(b),  9(b)  and  10(b).  With  the  number  of  cells 
increases  from  one  to  three,  the  proportion  of  heat  loss  from  walls 
for  the  whole  start-up  process  decreases  from  19.51%  to  7.04%,  and 


the  proportion  of  heat  absorbed  by  PEMFC  increases  from  78.35%  to 
90.58%.  The  results  indicate  that  increasing  the  number  of  cells  of 
stacks  improves  the  cold  start  performance. 

The  water  generation,  phase  change  and  transport  are  impor¬ 
tant  factors  affecting  the  cold  start  performance.  The  formed  water 
freezes  to  ice  during  a  cold  start  process  once  the  ionomer  is 
saturated.  Fig.  11  shows  the  water  content  evolution  in  diverse 
states  (ice  formation,  ionomer  absorption  and  evaporation)  of 
different  stacks  and  single  cell.  The  total  water  generation  rate  is 
the  same  for  the  individual  cells  of  different  stacks  and  single  cell  in 
the  galvanostatic  mode.  Therefore,  the  ice  formation,  ionomer 
absorption  and  evaporation  share  the  constant  water  generation. 
Among  them,  the  water  evaporation  (taken  away  by  gas)  is  insig¬ 
nificant  as  shown,  due  to  the  low  saturation  pressure  of  water  vapor 
at  subzero  temperatures.  The  discrepancies  between  the  different 
stacks  and  single  cell  are  subtle.  The  ionomer  in  the  stacks  with 
more  cells  can  absorb  more  water,  due  to  the  higher  saturation 
water  content  at  the  higher  temperature.  Fig.  11  explains  the  lower 
ice  formation  rate  of  the  stacks  with  more  cells,  and  as  the  number 
of  cells  further  increases,  it  is  suggested  that  this  effect  might  be 
more  notable. 

3.2.  Comparison  of  different  cells  in  three-cell  stack 

For  the  three-cell  stack,  the  individual  cells  are  defined  as  Cell  1, 
Cell  2  and  Cell  3,  representing  the  cell  on  the  cathode  side,  in  the 
middle  and  on  the  anode  side,  respectively,  as  demonstrated  in 
Fig.  1.  The  evolution  of  volume  averaged  cell  temperatures  and 
voltages  of  the  individual  cells  in  the  three-cell  stack  at  different 
current  densities  is  shown  in  Fig.  12.  It  can  be  noticed  that  no 
apparent  difference  exists  among  the  voltages  of  the  different 
individual  cells,  the  uniform  reactants  supply  is  one  of  the  reasons 
of  this  result.  The  evolution  of  temperature  for  different  cells  is 
similar  with  each  other,  with  slight  discrepancy  that  the  tempera¬ 
ture  of  Cell  2  is  the  highest  as  it  is  far  from  the  end  surfaces.  The 
temperature  of  Cell  3  is  higher  than  Cell  1  because  the  heat  source  is 
larger  on  the  anode  side  for  each  cell. 

3.3.  Comparison  of  the  single  cell  and  the  middle  cell  (Cell  2)  of 
three-cell  stack 

Since  Cell  2  (better  heat  insulation,  with  other  cells  on  its  both 
sides)  of  the  three-cell  stack  is  similar  to  the  individual  cells 
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Fig.  15.  Evolution  of  volume  averaged  cell  temperature  of  the  single  cell  and  Cell  2  of 
the  three-cell  stack  at  different  current  densities.  The  initial,  ambient  and  inlet  gas 

fraction  is  0. 
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between  the  two  end  cells  of  the  stacks  with  more  cells,  the  cold 
start  performance  of  Cell  2  is  representative  for  a  large  practical 
PEMFC  stack  with  many  individual  cells.  Therefore,  comparing  Cell 
2  of  the  three-cell  stack  and  the  single  cell  helps  understand  the 
difference  between  practical  stacks  and  single  cell.  Fig.  13  shows  the 
evolution  of  cell  voltages  and  ice  volume  fractions  in  the  cathode 
CLs  of  the  single  cell  and  Cell  2  of  the  three-cell  stack  at  different 
current  densities.  It  can  be  noticed  that  the  ice  formation  rate  of  the 
single  cell  is  higher  than  that  of  Cell  2.  The  cell  voltages  drop 
slightly  in  the  early  period  of  the  start-up  processes,  followed  by 
the  sharp  decrements  of  voltage  when  the  cathode  CLs  are  almost 
full  of  ice,  and  notable  difference  in  voltage  evolution  also  appears 
near  the  end  of  the  cold  start  processes.  For  the  single  cell,  since  the 
cathode  CL  is  blocked  by  ice  more  quickly,  the  sharp  voltage  drop 
happens  more  quickly  than  Cell  2  of  the  three-cell  stack. 

Fig.  15  shows  the  ice  volume  fraction  in  the  cathode  CLs  for  the 
single  cell  and  Cell  2  of  the  three-cell  stack  at  the  current  density  of 
0.1  A  cm  2.  Ice  is  first  formed  under  the  land  and  then  under  the  flow 
channel,  because  the  under  land  region  has  lower  temperature.  By 
comparing  Fig.  14(a)  and  (b),  it  can  be  found  that  for  the  single  cell, 
ice  is  formed  faster  than  Cell  2  of  the  three-cell  stack,  especially  in 
the  direction  from  the  inner  region  of  the  cathode  CL  to  BP.  This 
phenomenon  is  mainly  caused  by  the  different  temperatures. 

Fig.  15  shows  the  evolution  of  volume  averaged  cell  temperature 
of  the  single  cell  and  Cell  2  of  the  three-cell  stack  at  different 
current  densities.  It  can  be  noticed  that  the  average  temperature  of 
Cell  2  increases  faster  than  that  of  the  single  cell,  because  Cell  2  is 
thermally  insulated  by  the  outside  cells  while  the  single  cell  is 
exposed  to  ambient  directly.  More  detailed  feature  of  temperature 
distribution  can  be  found  in  Fig.  16.  This  figure  shows  part  of  the 
single  cell  and  Cell  2,  which  include  the  membrane,  cathode  CL, 
cathode  GDL,  cathode  flow  channel  and  BP,  the  scale  of  each 


component  is  adjusted  for  clear  recognition.  It  can  be  noticed  that 
the  temperature  distribution  of  Cell  2  is  more  uniform  than  the 
single  cell.  The  BP  of  the  single  cell  dissipates  heat  directly  to  the 
low  temperature  ambient,  and  because  the  thermal  conductivity  of 
BP  is  much  higher  than  the  other  components,  the  temperature 
distribution  in  the  BP  is  more  uniform,  and  large  temperature 
gradient  exists  in  the  MEA.  Furthermore,  with  the  high  thermal 
conductivity  of  the  BP,  in  the  MEAs,  the  temperature  is  lower  under 
the  land.  It  can  also  be  noticed  that  the  temperature  distribution  in 
the  MEA  is  more  uniform  in  Cell  2  than  the  single  cell,  because  the 
cooling  effect  of  the  BP  on  the  MEA  is  weaker  with  Cell  2.  It  can  also 
be  noticed  that  the  temperature  variation  is  insignificant  at  the 
membrane/CL  and  CL/GDL  interfaces,  due  to  the  similar  thermal 
conductivities  of  these  materials  (Table  1);  and  the  temperature 
variation  is  noticeable  at  the  GDL/BP,  GDL/channel  and  BP/channel 
interfaces,  due  to  the  different  thermal  conductivities  of  the 
materials  and  the  relatively  strong  convective  heat  transfer  in  the 
flow  channel. 

The  evolution  and  percentage  distribution  of  heat  generation 
and  loss  rates  of  the  single  cell  and  Cell  2  of  the  three-cell  stack 
through  various  ways  at  0.1  A  cm-2  are  shown  in  Fig.  17.  Comparing 
Figs.  17(a)  and  8(a)  show  that  the  heat  loss  rate  keeps  almost  zero 
within  the  first  3  s  for  Cell  2,  indicating  that  the  surrounding  does 
not  impact  the  inside  cells  of  stack  in  the  first  few  seconds.  The 
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Fig.  17.  Various  heat  generation  and  loss  rates  during  the  cold  start  process  of  Cell  2  of 
the  three-cell  stack  at  0.1  A  cm~2;  (a)  transient  variation;  (b)  percentage  distribution 
for  the  complete  cold  start  process.  The  initial,  ambient  and  inlet  gas  temperatures  are 
all  -20  °C,  the  initial  water  content  is  5,  and  the  initial  ice  volume  fraction  is  0. 


Fig.  18.  Evolution  of  volume  averaged  ice  volume  fraction  in  the  cathode  catalyst 
layers  of  the  individual  cells  in  the  three-cell  stack  at  different  current  densities.  The 
initial,  ambient  and  inlet  gas  temperatures  are  all  -20  °C,  the  initial  water  content  is  5, 
and  the  initial  ice  volume  fraction  is  0. 
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increment  of  heat  loss  from  walls  (the  virtual  middle  faces  of  BPs)  of 
Cell  2  is  much  lower  than  that  of  the  single  cell  due  to  the  better 
thermal  insulation  of  Cell  2.  Fig.  17(b)  shows  the  percentage 
distribution  of  various  heat  losses  integrated  though  the  whole 


starting  process.  The  heat  utilized  to  increase  the  temperature  is 
91.99%  of  the  total  heat  generation,  which  is  much  higher  than  the 
value  of  the  single  cell  of  76.67%.  The  better  heat  insulation  of  Cell  2 
indicates  that  most  of  the  cells  in  a  large  practical  PEMFC  stack  have 
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volume  fraction  during  the  start-up  processes  in  the  cathode  catalyst  layer  of  Cell  2  in  x-y  plane  in  the  center  of  the  stack  (25  mm  from  the 
A  cnr2.  The  initial,  ambient  and  inlet  gas  temperatures  are  all  -20  °C,  the  initial  water  content  is  5,  and  the  initial  ice  volume  fraction  is  0. 


Fig.  19.  Transient  distribution  of  ice  \ 
inlets)  at  the  current  density  of  0.05 , 
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better  cold  start  performance  than  a  single  cell,  while  the  cells 
located  outside  are  similar  to  a  single  cell  due  to  the  direct  heat  loss 
to  the  ambient.  Therefore,  strategies  to  assist  cold  start  should  focus 
on  the  outside  cells. 

3.4.  Effect  of  current  density 

Fig.  18  shows  the  evolution  of  volume  averaged  ice  volume 
fraction  in  cathode  CL  of  each  cell  of  the  three-cell  stack  at  different 
current  densities,  and  no  noticeable  difference  can  be  found 
between  the  different  cells  at  the  same  current  density.  However, 
the  influence  of  current  density  is  significant.  A  higher  current 
density  results  in  a  higher  ice  formation  rate,  and  as  Fig.  18(b) 
shows,  at  the  current  density  of  0.1  A  cm-2,  ice  starts  to  be  formed 
at  about  0.8  s,  while  at  0.05  A  cm-2,  ice  appears  at  around  2  s.  The 
reason  is  that  with  a  lower  water  generation  rate  (lower  current 
density),  the  ionomer  has  more  time  to  absorb  the  product  water, 
on  the  other  hand,  at  a  higher  current  density,  the  ionomer  in  the 
cathode  CL  gets  saturated  more  quickly. 

Fig.  19  shows  the  transient  ice  volume  fraction  distribution  in 
the  cathode  CL  in  Cell  2  of  the  three-cell  stack  at  0.05  A  cm-2,  as 
a  comparison,  Fig.  14(b)  shows  the  ice  distribution  at  the  same 
location  while  at  a  higher  current  density  (0.1  A  cm-2).  It  can  be 
found  that  in  both  cases  the  ice  appears  first  under  the  land. 
However,  at  the  lower  current  density,  ice  is  formed  faster  at 
a  section  close  to  the  BP  side  than  the  membrane  side;  on  the  other 
hand,  at  the  higher  current  density,  ice  is  formed  at  a  section  close 
to  the  membrane  side.  The  different  ice  distributions  at  different 
current  densities  are  caused  by  start-up  duration.  With  a  lower 
current  density,  the  membrane  water  absorption  effect  is  more 
significant  due  to  the  longer  duration,  therefore,  more  water  at  the 
sections  close  to  the  membrane  can  be  absorbed  by  the  membrane; 
on  the  other  hand,  such  effect  is  less  significant  at  a  higher  current 
density,  and  the  sections  close  to  the  membrane  side  form  ice  first 
because  of  the  higher  reaction  rate. 

4.  Conclusion 

To  understand  the  details  of  proton  exchange  membrane  fuel 
cell  (PEMFC)  stack  cold  start  processes,  a  three-dimensional 
multiphase  stack  model  has  been  developed  and  numerical  simu¬ 
lations  have  been  conducted.  It  is  found  that  notable  voltage 
discrepancies  between  different  stacks  and  single  cell  occur  when 
the  ice  almost  fully  blocks  the  cathode  catalyst  layer  (CL),  and  for 
the  stacks  with  more  cells,  the  voltage  decreases  more  slowly  due 
to  the  lower  ice  formation  rate.  The  temperature  increases  faster  for 
a  stack  with  more  cells,  and  a  higher  temperature  can  be  reached  at 
the  end  of  the  cold  start  process.  No  apparent  difference  in  voltage 
exists  among  the  different  individual  cells  in  a  stack  when  the 
reactant  gases  are  evenly  supplied  to  each  cell.  The  temperature  of 
the  individual  cell  in  the  middle  of  a  stack  is  higher  and  more 
evenly  distributed  than  those  on  the  sides  and  single  cells,  due  to 
weakened  cooling  effect  of  the  bi-polar  plate  (BP)  on  the  membrane 
electrode  assembly  (MEA),  and  the  ice  formation  rate  is  also  lower 
in  the  middle  cell.  The  temperature  of  the  individual  cell  in  a  stack 
on  the  anode  side  is  higher  than  on  the  cathode  side  due  to  the 
higher  ohmic  heat  in  the  anode.  At  a  lower  current  density,  the  ice 
in  the  cathode  CL  is  formed  faster  at  the  section  close  to  the  BP,  and 
it  is  close  to  the  membrane  at  a  higher  current  density. 
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